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We present an experimental comparison between different iterative ghost imaging algo- 
rithms. Our experimental setup utilizes a spatial light modulator for generating known 
random light fields to illuminate a partially-transmissive object. We adapt the weighting 
factor used in the traditional ghost imaging algorithm to account for changes in the 
efficiency of the generated light field. We show that our normalized weighting algorithm can 
match the performance of differential ghost imaging. © 2012 Optical Society of America 

OCIS codes: (030.4280) Noise in imaging systems; (030.6140) Speckle; (110.1650) Coherence imag- 
ing; (200.1130) Algebraic optical processing. 

References 

1. R. S. Bennink, S. J. Bentley, and R. W. Boyd, " "Two-photon" coincidence imaging with a classical 
source," Phys. Rev. Lett. 89, 113601 (2002). 

2. A. Gatti, E. Brambilla, M. Bache, and L. A. Lugiato, "Ghost imaging with thermal light: Comparing 
entanglement and classical correlation," Phys. Rev. Lett. 93, 093602 (2004). 

3. A. Gatti, E. Brambilla, M. Bache, and L. A. Lugiato, "Correlated imaging, quantum and classical," 
Phys. Rev. A 70, 013802 (2004). 

4. A. Valencia, G. Scarcelli, M. D'Angelo, and Y. Shih, "Two-photon imaging with thermal light," Phys. 
Rev. Lett. 94, 063601 (2005). 

5. F. Ferri, D. Magatti, A. Gatti, M. Bache, E. Brambilla, and L. A. Lugiato, "High-resolution ghost image 
and ghost diffraction experiments with thermal light," Phys. Rev. Lett. 94, 183602 (2005). 

6. J. H. Shapiro, "Computational ghost imaging," Phys. Rev. A 78, 061802 (2008). 

7. Y. Bromberg, O. Katz, and Y. Silberberg, "Ghost imaging with a single detector," Phys. Rev. A 79, 
053840 (2009). 

8. M. Duarte, M. Davenport, D. Takhar, J. Laska, T. Sun, K. Kelly, and R. Baraniuk, "Single-pixel imaging 
via compressive sampling," IEEE Signal Processing Magazine 25, 83-91 (2008). 

9. F. Ferri, D. Magatti, L. A. Lugiato, and A. Gatti, "Differential ghost imaging," Phys. Rev. Lett. 104, 
253603 (2010). 

10. J. Goodman, Statistical Optics (Wiley, 2000). 

11. S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University Press, 2004). 

12. D. L. Donoho and Y. Tsaig, "Fast solution of li-norm minimization problems when the solution may be 
sparse," IEEE Transactions on Information Theory 54, 4789-4812 (2008). 



1 



13. O. Katz, Y. Bromberg, and Y. Silberberg, "Compressive ghost imaging," Appl. Phys. Lett. 95(13), 
131110 (2009). 



1. Introduction 

Classical ghost imaging (GI) [I|j5] uses a series of random light patterns to illuminate an unknown object. 
For each pattern the reflected or transmitted light is measured using a single element detector. The series 
of single element measurements, combined with the known light patterns is used to deduce the object. In 
some systems the random light pattern is produced as a time varying laser speckle, and a beam splitter is 
used to illuminate both the unknown object and a reference camera, with which the pattern is recorded. 
Subsequently, the need for the beam splitter and camera has been removed by implementing a spatial light 
modulator (SLM) to produce a random, but known, pattern thereby reducing the number of components in 
the system necessary for GI experiments pju]. This latter approach is known as computational GI and in 
terms of the experimental arrangement is closely related to the field of single pixel cameras [§]. 

In all approaches to GI an algorithm is employed to deduce the object using the series of measurements 
from the single element detector and either the recorded or computationally predicted random patterns. The 
algorithms employed fall into two categories, iterative ones that give a refined estimate of the object after 
every new light pattern and measurement, and inversion ones which infer an object based on the entire series 
of patterns and measurements. 

Iterative algorithms use the measured signal to derive a weighting factor to the corresponding pattern that 
is then added to the iterative estimate of the object. In this paper we compare a number of these iterative 
algorithms within the context of computational GI. The algorithms we consider are traditional GI (TGI) and 
differential GI (DGI) (9J. In a computational GI setup, TGI uses a weighting factor equal to the signal from 
the detector whereas DGI utilizes a weighting factor that depends on fluctuations in the measured signal 
and uses an additional detector to give a normalization. Beyond these two algorithms we introduce a variant 
of the TGI algorithm, normalized GI (NGI), which we show can match the performance of DGI. 

Key to all these algorithms is that the changes in the measured signal should arise from the overlap of 
the known random pattern with the unknown object. Obviously other sources of signal change are possible; 
including fluctuations arising from changes in the source intensity and changes in the efficiency with which 
the pattern is imprinted. These later sources of noise scale with the signal level and hence become more 
significant when the signal is high. 

2. Experimental setup 

The experimental setup is shown in Figure [I] Here a random light pattern is generated from a simulated 
superposition of plane waves using random numbers, which is then sent to an SLM to produce a synthesized 
speckle field. The SLM has 512x512 pixels in the window of size 3.584 x 3.584 mm. We pass a collimated laser 
of wavelength A = 632.8 nm through a polarizing beam splitter and a half-wave plate, before illuminating 
the SLM window. The speckle field is generated by modulation of the SLM and the returning light field is 
then magnified by a simple telescope system consisting of 150 mm and 450 mm biconvex lenses. The object 
is located at the focus plane of the 450 mm lens, which is also the image plane of the SLM window. A 50 : 50 
beam splitter is placed before the object in order to split the speckle field into two beams; the object beam 
(I(xs)) and the reference beam (I(xr)). The object beam illuminates the object and is then collected by a 
bucket detector, thus providing an computational GI setup. The additional reference beam for monitoring 
the light differentiates our system from previous experimental computational GI configurations. Since we are 
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generating a computer hologram that is then sent to the SLM to create the speckle field, we can therefore 
predict the light field at the reference arm, negating the demand for a CCD camera, and requiring only a 
second bucket detector. It should be noted that for TGI based on our computational GI setup, only the 
object bucket detector is needed. The additional bucket detector in the reference arm is only required for 
NGI and DGI. Light intensities detected by the object and reference bucket detectors are indicated by S 
and R respectively, and the speckle field is described by I(x,y). As we use a 50 : 50 beam splitter, it is 
understood that I(x,y) = 2I(x s ,ys) = 2I(x R ,y R ). 




150 mm 450 mm 



Fig. 1. Computational ghost imaging setup used in the experiment. A spatial light modulator 
(SLM) is used to generate a random speckle field, as described in the text and a beam splitter 
(BS) is used to measure a reference signal R on a bucket detector before the object. The 
signal, S, is measured on a bucket detector which collects the light transmitted after the 
object. 



3. Iterative ghost imaging algorithms 

In all iterative GI techniques, the transmitting object located after the beam splitter, BS2, is reconstructed 
by correlating the speckle field intensity measured at S and R, then adding together each successive frame 
with a suitable weighting factor. The transmitted light power detected after the object can be expressed as 

S= I(xs,ys)T(x s ,ys)dx s dys, (1) 

where the laser area is A\ and T(xs, ys) is the (intensity) object transmission function, while the background 
reference is expressed as 

R = J I(x R ,y R )dx R dy R . (2) 

3. A. Traditional Ghost Imaging 

In TGI, the reconstruction result of the object, 0(x,y) is retrieved from the correlation between S and 
I(x,y). We define for each iteration, z, the contribution to the reconstruction to be 7 

Oi(x, y) = (S- (S)) (I(x, y) - (I(x, y))) , (3) 
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where < . >= denotes an ensemble average for M iterations. We obtain the final reconstruction by 

averaging over all iterations such that 0(x,y) = (Oi(x,y)). It is easy to understand the reconstruction as 
being derived from the weighted sum of the speckle field for each measurement. Therefore S is the weight 
for the speckle field for each measurement. One drawback of using this algorithm is that the reconstruction 
is heavily weighted to the size of the signal S and is thus susceptible to fluctuations in the generated light 
field. These fluctuations can arise from either changes to the laser power or the efficiency of the SLM in 
computational GI. 

3.B. Differential ghost imaging 

Differential GI 9 , first performed by Ferri et al, utilizes a second bucket detector to extract a reference signal 
which is used in the reconstruction to weight the speckle field based on the average transmission signal relative 
to the average reference signal. Similarly, each contribution to the reconstruction can be expressed as 

Oi(x, y)=(s- (/(*, y) - (I(x, y))) . (4) 

Thus we obtain the final result by summing for all iterations. We observe the second term in brackets on the 
right hand side of Eq.[3]and Eq.|4]are both identical however the first term in brackets of Eq.[4]is now weighted 
according to the average value of 5, which is normalized to the average value of R. As demonstrated in [9] the 
DGI algorithm improves by order of magnitude the SNR of the measurement with respect to TGI. Moreover, 
a key difference from TGI, it is no longer sensitive to other sources of noise. For example, fluctuations in the 
laser power or changes to the SLM efficiency will affect both the reference signal and the transmitted signal, 
and thus the contribution to the reconstruction will be weighted more appropriately. 

3.C. Normalized ghost imaging 

3.C.I. Normalized ghost imaging with two detectors 

As seen in Eq. [4j larger values of S measured by the bucket detector results in a greater weight for that 
particular speckle field, therefore external noise sources can still affect the overall reconstruction. There exists 
another iterative algorithm which instead normalizes each individual measurement 5, as well as the running 
average, according to the reference signal R, resulting in an arguably more intuitive approach for dealing 
with time varying noise sources. We call this approach normailized GI (NGI). The algorithm used to describe 
each contribution to the reconstruction in NGI is given by 

Oi(x, 2/) = (| " H) {Ox, y) - (I(x, y))) , (5) 

where we have assumed ~ for a large number of measurements. By considering equations jij and J^J 
we can summarize the difference between the two algorithms as 

(0(x, y) NGI ) = -^y (0(x, y) D Gi) • (6) 

3.C.2. Normalized ghost imaging with a single detector 

In a computational GI setup, we can show that the additional detector used to measure the reference signal in 
DGI and NGI can instead be estimated based on the known light field reflected from the SLM and the average 
measured signal S for an arbitrary number of previous iterations. Calculating R negates the requirement 
for an additional detector, whilst improving the performance of the reconstruction compared to TGI, thus 
single-detector NGI (SNGI) is identical to the TGI experimental setup, with only a modified algorithm. 
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3.D. Signal-to-noise ratio analysis 

To make a quantitative comparison between the NGI and the existing algorithms, we adopt a similar approach 
as used by Ferri et al and investigate the theoretical contribution to the signal-to- noise ratio (SNR) for objects 
with varying transmission functions. In 9 the authors express the average quantity of Eq.[4]in terms of the 
object transmission fluctuation 5T(x, y) = T(x, y) — T, 

(0(x,y) DGI ) = A~ s (I) 2 ST(x,y), (7) 

where A s is the average speckle area and T = j Ai (I(x,y))T(x,y)dxdy/ J A (I(x,y)) dxdy is the average 
transmission function of the object. Note that Eq.[7|is obtained under the assumptions of uniform illumination 
(the average speckle beams are constant over their area) and perfect resolution (the speckle area is much 
smaller compared to features of the object). The corresponding signal of DGI can be defined as 

(A(0 DGI )) 2 = A; 2 (I) 4 (AT) 2 , (8) 

where AT is the variation of the object transmission function to be detected. Similarly, using Eq. [6j we can 
express the signal of NGI as 

(a(o^ g/ )) 2 = a; 2 ^(at) 2 . (9) 

{H) 

The speckle patterns used in our experiment exhibit complex-Gaussian behaviour, such that the intensity 
is exponentially distributed (see Fig. [2]), and the noise associated to the measurement of 0(x,y) can be 
expressed as 

(S0 2 (x, y)) = (0(x, yf) - (0(x, y)) 2 , (10) 

for which it can be shown that (0(x, y)) =0, thus the second term on the right hand side (RHS) of in Eq.[lQ| 
may be omitted. Again, under the assumptions of uniform illumination and perfect resolution, the noise of 
DGI can be expressed as 

(Olc^vAsAtil) 4 ^, (11) 
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where ST 2 — T 2 — T and T 2 = J A (I(x,y))T 2 (x,y)dxdy/ J A (I(x,y)) dxdy. Using linearization we can 
write 

S„W( 1+ M-S*\ (12) 
where SS and SR are the zero-mean deviation of S and R, thus the noise of NGI is shown to be 

(O^^AsA^ST 2 -. (13) 
{H) 

Finally, we show that the SNR contribution for NGI is 

M AT 2 

SNR NGI = SNRjjgi = „ ( 14 ) 

-lv speckle 01 

where N s = A\jA s is the number of speckles in the field. The SNR contribution for NGI is found to be 
identical to that of the DGI algorithm derived in 9 . For comparison the SNR contribution for TGI was 
shown to be 

M AT 2 

S ™™' = N^TT- (15) 

Therefore we can examine the difference between the NGI (or DGI) and TGI algorithms by obtaining the 
ratio of SNR calculations, given as 

2 

SNR NGI _ T 
SNRtgi T2-T 2 ' 

As highlighted by Ferri et al, the difference is always greater than 1 and dependent only upon the variation 
in the object transmission function. 

4. Experiment results 

We generated a series of random speckle patterns using an SLM by simulating the interference of many plane 
waves on a computer. The real and imaginary amplitude components and the wave vector k of each simulated 
plane wave is Gaussian distributed. Figure[2] shows a typical example of the speckle patterns generated on the 
SLM and the exponentially distributed intensity for many patterns, implying that the speckle hologram has 
complex- Gaussian statistics, thereby a good approximation for real speckle fields [To]. A binary transmissive 
object, 5mm x 5mm in size, is located after a 3x magnification telescope in the image plane of the SLM. 
Since we know both the object and the random speckle field projected to the SLM, we are able to simulate 
the expected results for comparison with our experiment. Experimental and simulated reconstruction results 
after 10000 iterations are shown in Fig. [3J The simulated reconstruction is produced assuming no external 
noise sources. The partially transmissive object used is indicated in the bottom right of Fig. [3j It is clear 
that the DGI and NGI algorithms provide very similar results, as predicted from the theory, and both show 
improved background subtraction compared to TGI. 

Compared with the traditional computational GI setup, the NGI algorithm requires a reference bucket 
detector. However, as discussed in section |3.C.2[ the advantage of computational GI means that we can 
replace this bucket detector with a virtual reference detector generating a simulated R. Thus we can negate 
the requirement for the reference detector and return the system to a true single element camera, which we 
call single-detector NGI (SNGI). The two major factors that dominate the value of R are from the different 
speckle patterns displayed on the SLM and fluctuations of the incident laser power. We can computationally 
predict changes to the value of R due to the speckle pattern, whereas fluctuations of the laser power can be 
simulated by using a rolling average for a particular series of S measurements. The bottom row in Fig. [3] 
shows the experimental results for reconstructing the object using the SNGI algorithm. We observe similar 



6 



algorithm 


experiment 


simulation 






























TGI 




























m 
































DGI 
















* 












(2 




NGI 


IIP® 

'•" • ■*■" i S ! n 






IPI 






Ha 












(3 




SNGI 


























































iterations 


10 


100 


1000 


10000 


object 



Fig. 3. Experimental results (middle column) for TGI, DGI and NGI reconstruction algo- 
rithms as they evolve (10, 100, 1000 and 10000 iterations from left to right, respectively) 
with the corresponding simulated results (right column). The transmissive object is shown 
in the lower right. The bottom row shows the evolution for reconstructing the object with 
the NGI algorithm using a single detector and predicting the reference signal i?, termed 
here the SNGI algorithm. 

results compared with DGI and NGI algorithms indicating an improved performance compared with the 
TGI algorithm for single element camera. 

To demonstrate the effect of object transmission function on the performance of NGI compared with TGI 
and DGI algorithms we used a similar experimental approach to that in Ref. 9 . By scanning a knife edge 
(located in the image plane of the SLM, as before) across the speckle field in well defined steps (for which 
AT = 1), we measured the SNR's for the final object reconstruction obtained after 5000 random speckle 
iterations. The beam size used was 10x10 mm and the speckle size at the plane of the object was found to 
be S s ~ 90 /im, providing around N s ~ 12500 speckles. The experimental results and theoretical predictions 
for the SNR's of each iterative algorithm are shown in Fig. [4] Note that the y-axis has been normalized to the 
number of iterations. We observe close quantitative agreement between the theory and the measurements. 
The results indicate that for low transmissive objects, all algorithms reconstruct with similar SNR, while for 
more transmissive objects the DGI and NGI algorithms become more efficient in comparison to TGI due 
to the differential nature of the reconstruction. Furthermore, we observe that when using a single detector, 
SNGI is a more efficient algorithm for reconstructing objects of all transmissions compared to TGI. We 
observe that for increasing transmissive objects SNGI becomes less efficient than NGI, for which the reason 
is the subject of ongoing research. Similar to 9 , we find a systematic discrepancy between the experimental 
results of TGI and the theoretical predictions. 
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Fig. 4. Signal-to-noise ratio's for DGI, NGl, SNGI and TGI versus transmitting area. Trans- 
mitting ratio is defined as the ratio between the transmitting area of the object and the 
area of the speckle field. 



5. Normalization in matrix inverse algorithms 

5. A. Introduction to matrix inverse algorithms and compressive sensing 

As an alternative to the iterative techniques discussed above, we can choose to record all the signals for a 
complete set of speckle patterns and then treat the image reconstruction as one of matrix inversion. The 
series of M speckle patterns, each containing N pixels can be represented by a M x N matrix. If the object 
is also represented as an N element column vector, then the vector containing the measured signals is a M 
element vector. This relationship is expressed as 



(17) 



In the case where the number of speckle patterns equals the number of pixels then the M x N matrix is 
square, such that its inverse can be calculated and the object vector determined. However when M < N and 
or TV is large, the system is ill-conditioned and calculating the inverse of the matrix is not straightforward. 
Problems of this type are wide spread in physics and techniques for solving them have been developed. Within 
our system the appeal is to reconstruct the image of N pixels from M measurements where M < N . That 
this is possible is based on the fact that natural images are sparse and the reconstruction can be obtained 
by solving a convex optimization problem |11 , which is a generalization of a linear least squares problem. 
In contrast to iterative methods, compressive GI (CGI) needs to take all measurements, represented here, in 
some compressible basis (in this case a discrete cosine transform which has been applied to each row of the 
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M x N matrix). Solving the convex optimization problem requires minimizing the £1 norm 12 . 
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5.B. Normalized compressive ghost imaging 

By normalizing the measured object signal relative to the reference signal as performed above, such that 



S' = S/R, we can apply the CGI technique 13 to reconstruct our object. Equation 
for normalized CGI (NCGI) as 
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can then be written 



(18) 



Performing both NCGI and CGI analyses using the same experimental data (acquired using the exper- 
imental setup in FigjlJ we obtain the reconstruction in Fig. [5] We observe a clear improvement using the 
NCGI algorithm compared to the CGI algorithm, manifest as an increased SNR value. The efficiency with 
which NCGI can reconstruct sparse images over CGI is determined by the level of noise in the system. We 
find that when there is no system noise present, both reconstructions are essentially identical. Thus the main 
improvement in employing NCGI over CGI with the additional reference detector is the ability to protect 
the reconstruction from time varying noise sources. 

6. Conclusion 

In conclusion we have compared different iterative GI methods to reconstruct an object and studied a 
new GI algorithm, which we call normalized GI (NGI). The performance of the differential GI (DGI) and 
NGI algorithms show good quantitative agreement as predicted by the theoretical foundations that support 
them. Our results indicate that by normalizing the measured signal relative to a reference signal, a more 
appropriate weighting factor is applied to the ensemble average of the estimated object, compared to the 
traditional GI (TGI) algorithm. Our analysis of the measured SNR and the object transmission shows a 
significant improvement for more transmissive objects in comparison to TGI. Furthermore, we have shown 
it is possible to apply normalization to systems with a single detector, SNGI, by estimating the reference 
signal. We have also investigated normalization within a compressive matrix inversion method, showing 
similar results to an non-normalized algorithm but with enhanced noise suppression. We believe the NGI 
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algorithm will be a useful resource for imaging where alternative techniques are required in the future. 
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